# ------------------------------------------------------------------------------#
# Description: Create Table 3: Peer effects on adoption of any RainWise component (cistern, rain garden or both)
# Project:     Peer Effects in Voluntary Environmental Policies:
#              An Application to Urban Water Quality
# Author:      Daniel A. Brent | dab320@psu.edu
# Created:     2025-06-11
# ------------------------------------------------------------------------------#

# Clear environment ------------------------------------------------------------
rm(list = ls())
gc()

# Load required packages -------------------------------------------------------
library(pacman)
p_load(data.table, fixest, modelsummary, stringr)

options(scipen = 999)

# Load adoption data -----------------------------------------------------------
load("data/adoptions/adoptions.RData")

# Merge peer adoption buckets --------------------------------------------------
load("data/bucket/peer_adopt_buckets_any.RData")
dat <- merge(adopt, peer.adopt.bucket.any, by = c("pin", "year"))
rm(adopt, peer.adopt.bucket.any)

load("data/bucket/peer_adopt_buckets_rg.RData")
dat <- merge(dat, peer.adopt.bucket.rg, by = c("pin", "year"))
rm(peer.adopt.bucket.rg)

load("data/bucket/peer_adopt_buckets_any_rg.RData")
dat <- merge(dat, peer.adopt.bucket.any.rg, by = c("pin", "year"))
rm(peer.adopt.bucket.any.rg)

load("data/bucket/peer_adopt_buckets_cs.RData")
dat <- merge(dat, peer.adopt.bucket.cs, by = c("pin", "year"))
rm(peer.adopt.bucket.cs)

load("data/bucket/peer_adopt_buckets_both.RData")
dat <- merge(dat, peer.adopt.bucket.both, by = c("pin", "year"))
rm(peer.adopt.bucket.both)

# Merge peer eligibility buckets -----------------------------------------------
load("data/bucket/peer_elig_buckets_cs.RData")
dat <- merge(dat, peer.buckets.cs, by = c("pin", "year"))
rm(peer.buckets.cs)

load("data/bucket/peer_elig_buckets_rg.RData")
dat <- merge(dat, peer.buckets.rg, by = c("pin", "year"))
rm(peer.buckets.rg)

# Merge total peers ------------------------------------------------------------
load("data/bucket/total_peers.RData")
dat <- merge(dat, total.peers, by = "pin")
rm(total.peers)

# Merge parcel data ------------------------------------------------------------
load("data/parcel_build/parcel_build.RData")
dat <- merge(dat, parcel.build[, .(pin, zip5, sub.area, sqft, lot, beds, baths,
                                   yearbuilt, val.land, val.improve,
                                   ren.pre90, ren.90.99, ren.00.09, ren.10,
                                   water.prob)], by = "pin")
rm(parcel.build)

# Merge census data ------------------------------------------------------------
load("data/census/all_parcels_census_2010.RData")
parcels.2010[, blk := as.numeric(paste0(block, blkgrp))]
dat <- merge(dat, parcels.2010[, .(pin, tract, blkgrp, block, blk)])
rm(parcels.2010)

# Scale peer adoption variables ------------------------------------------------
dat[, peer.adopt.any.100.scale     := peer.adopt.any.100 / 100]
dat[, peer.adopt.any.rg.100.scale  := peer.adopt.any.rg.100 / 100]
dat[, peer.adopt.rg.100.scale      := peer.adopt.rg.100 / 100]
dat[, peer.adopt.cs.100.scale      := peer.adopt.cs.100 / 100]
dat[, peer.adopt.both.100.scale    := peer.adopt.both.100 / 100]

# Keep post-2010 only ----------------------------------------------------------
dat <- dat[year > 2010]

# Create average eligibility variable ------------------------------------------
dat[, elig.peers.avg.100 := elig.peers.cs.100 + elig.peers.rg.100]

# Table 3: Main Results --------------------------------------------------------------------------

## OLS --------------------------------------------------------------------------
dat[, peer.adopt.any.100.ols := peer.adopt.any.100]

m1 <- feols(adopt.any * 100 ~ peers.100 + peer.adopt.any.100.ols | 
              year + blkgrp,
            cluster = "blkgrp",
            data = dat[elig.cs == 1 & post.rainwise == 0])

## First Stage ------------------------------------------------------------------
m2 <- feols(peer.adopt.any.100 ~ peers.100 + elig.peers.cs.100 + elig.peers.rg.100 | 
              year + blkgrp,
            cluster = "blkgrp",
            data = dat[elig.cs == 1 & post.rainwise == 0])

## IV - Multi -------------------------------------------------------------------
m3 <- feols(adopt.any * 100 ~ peers.100 | 
              year + blkgrp | 
              peer.adopt.any.100 ~ elig.peers.cs.100 + elig.peers.rg.100,
            cluster = "blkgrp",
            data = dat[elig.cs == 1 & post.rainwise == 0])

# IV - Single ------------------------------------------------------------------
m4 <- feols(adopt.any * 100 ~ peers.100 | 
              year + blkgrp | 
              peer.adopt.any.100 ~ elig.peers.avg.100,
            cluster = "blkgrp",
            data = dat[elig.cs == 1 & post.rainwise == 0])

## Diagnostics ------------------------------------------------------------------
fitstat(m3, ~ivf + ivwald + ivfall + ivf1 + wh + sargan)

varnames <- c(
  'peer.adopt.any.100.ols' = 'Peer Adoptions',
  'elig.peers.cs.100' = '\\# Eligible Peers (CS)',
  'elig.peers.rg.100' = '\\# Eligible Peers (RG)',
  'peer.adopt.any.100' = '$\\widehat{\\text{Peer Adoptions}}$',
  'peer.adopt.any.100.scale' = '$\\widehat{\\text{Peer Adoptions}}$'
)

# Load IV diagnostics ----------------------------------------------------------
load("data/iv/iv_diag_table3_4.RData")
names(iv_diag_table3_4) <- c("table3.iv.mult", "table3.iv", "table4.any.rg",
                             "table4.cs", "table4.rg", "table4.both")

eff.f.mult   <- round(iv_diag_table3_4[["table3.iv.mult"]]$F_stat[4], 2)
eff.f.sing   <- round(iv_diag_table3_4[["table3.iv"]]$F_stat[4], 2)
ar.ci.mult   <- paste0("[", round(iv_diag_table3_4[["table3.iv.mult"]]$AR$ci[1], 2), ", ",
                       round(iv_diag_table3_4[["table3.iv"]]$AR$ci[2], 2), "]")
ar.ci.sing   <- paste0("[", round(iv_diag_table3_4[["table3.iv"]]$AR$ci[1], 2), ", ",
                       round(iv_diag_table3_4[["table3.iv"]]$AR$ci[2], 2), "]")
tf.ci.sing   <- paste0("[", round(iv_diag_table3_4[["table3.iv"]]$tF[6], 2), ", ",
                       round(iv_diag_table3_4[["table3.iv"]]$tF[7], 2), "]")
tf.cf.sing   <- round(iv_diag_table3_4[["table3.iv"]]$tF[2], 2)

## Preview etable ---------------------------------------------------------------
etable(m1, m2, m3, m4,
       fixef_sizes = TRUE, depvar = FALSE, se.below = TRUE, dict = varnames,
       drop = "peers.100",
       fitstat = c("n", "sargan", "wh"),
       extralines = list(
         "_F"       = c("", "", eff.f.mult, eff.f.sing),
         "_AR CI"   = c("", "", ar.ci.mult, ar.ci.sing),
         "_tF cF"   = c("", "", "", tf.cf.sing),
         "_tF CI"   = c("", "", "", tf.ci.sing)
       ))

## Format LaTeX style -----------------------------------------------------------
style_tex <- style.tex(
  tablefoot.value = '',
  model.title = '',
  depvar.title = '',
  var.title = '\\midrule',
  stats.title = "\\midrule"
)

setFixest_etable(
  drop = "peers.100",
  fitstat = c("n"),
  digits = 3, digits.stats = 3,
  depvar = FALSE,
  style.tex = style_tex
)

## Export Table 4 to LaTeX ------------------------------------------------------
etable(m1, m2, m3, m4, tex = TRUE,
       fixef_sizes = TRUE, depvar = FALSE, se.below = TRUE, dict = varnames,
       signif.code = c("*" = 0.1, "**" = 0.05, "***" = 0.01),
       fitstat = c("n", "sargan", "wh"),
       headers = list("_ " = c("OLS", "First Stage", "IV", "IV")),
       extralines = list(
         "_F"       = c("", "", eff.f.mult, eff.f.sing),
         "_AR CI"   = c("", "", ar.ci.mult, ar.ci.sing),
         "_tF cF"   = c("", "", "", tf.cf.sing),
         "_tF CI"   = c("", "", "", tf.ci.sing)
       ),
       file = "output/tables/table_3.tex", replace = TRUE)
